from EOS_collection import antoine, G_ex, gamma, NRTL
import numpy as np
from scipy.optimize import minimize
import math
import matplotlib.pyplot as plt

#2024.08.05
#NRTL方程的拟合
#2024.12.11
#tau12,tau21与温度简单关联
#
#
#基础参数以甲醇和水为例
T1 = np.array([365.54, 364.88, 364.12, 363.19, 362.13, 361.21, 360.05, 359.04, 358.02, \
               357.53, 356.70, 355.56, 354.35, 353.91, 353.24, 352.32, 351.488, 350.76, \
               350.02, 348.61, 348.00, 347.21, 345.98, 344.33, 342.85, 340.56])
x_1 = np.array([0.0484, 0.0539, 0.0607, 0.0696, 0.0806, 0.0911, 0.1056, 0.1197, \
                0.1354, 0.1436, 0.1586, 0.1815, 0.2092, 0.2202, 0.2380, 0.2648, \
                0.2918, 0.3166, 0.3442, 0.4015, 0.4282, 0.4644, 0.5240, 0.6091, 0.6899, 0.8213])
y_1 = np.array([0.2704, 0.3009, 0.3213, 0.3492, 0.3866, 0.4126, 0.4492, 0.4762, \
                0.4885, 0.5047, 0.5323, 0.5677, 0.5969, 0.6031, 0.6237, 0.6519, \
                0.6681, 0.6923, 0.7197, 0.7431, 0.7443, 0.7613, 0.7860, 0.8231, 0.8633, 0.9221])
A1 = 7.19736
B1 = 1574.99
C1 = -34.29
A2 = 7.07404
B2 = 1657.46
C2 = -46.13
atm = 101.325
x_1_exp = x_1
y_1_exp = y_1
x_2 = 1 - x_1
y_2 = 1 - y_1

#压力，气体常数
Pa = np.ones_like(T1) * atm
R = 8.314
alpha = 0.3

#控制输出的精度
np.set_printoptions(precision=3)

#计算饱和蒸汽压
Pvs_1 = antoine (A1, B1, C1, T1)
print(f"组分1的饱和蒸汽压{Pvs_1}")
Pvs_2 = antoine (A2, B2, C2, T1)
print(f"组分2的饱和蒸汽压{Pvs_2}")

#计算活度系数
gamma_1 = gamma(x_1, y_1, Pa, Pvs_1)
gamma_2 = gamma(x_2, y_2, Pa, Pvs_2)
print(f"组分1的活度系数{gamma_1}")
print(f"组分2的活度系数{gamma_2}")

#计算过剩吉布斯自由能
G_excess = G_ex (x_1, x_2, gamma_1, gamma_2, T1)
# print(f"过剩吉布斯自由能{G_excess} J/mol")

#拟合NRTL方程
def NRTL_fit(a, alpha):
    tau12 = a[0] + a[1] / T1
    tau21 = a[2] + a[3] / T1
    G21 = np.exp(-alpha * tau21)
    G12 = np.exp(-alpha * tau12)
    gamma_1_cal = NRTL(x_1, x_2, tau12, tau21, alpha)
    gamma_2_cal = NRTL(x_2, x_1, tau21, tau12, alpha)
    f = (gamma_1_cal - gamma_1) ** 2 + (gamma_2_cal - gamma_2) ** 2
    return np.sum(f) 

#拟合NRTL方程
guess = np.ones(4)
a = minimize(NRTL_fit, guess, alpha, method = 'L-BFGS-B').x
tau12 = a[0] + a[1] / T1
tau21 = a[2] + a[3] / T1
# print(f"a={a}")
# print(f"tau12 = {tau12}")
# print(f"tau21 = {tau21}")

#验算活度系数
gamma_1 = NRTL(x_1, x_2, tau12, tau21, alpha)
gamma_2 = NRTL(x_2, x_1, tau21, tau12, alpha)
# print(f"验算组分1的活度系数{gamma_1}")
# print(f"验算组分2的活度系数{gamma_2}")

#计算相图
x_1 = np.linspace(0, 1, 20)
x_2 = 1 - x_1
y_1 = np.ones_like(x_1)
T2 = np.ones_like(x_1) * 300
Pa = np.ones_like(x_1) * atm

#求温度
def equilibrium(T2):
    tau12 = a[0] + a[1] / T2
    tau21 = a[2] + a[3] / T2
    gamma_1 = NRTL(x_1, x_2, tau12, tau21, alpha)
    gamma_2 = NRTL(x_2, x_1, tau21, tau12, alpha)
    Pvs_1 = antoine (A1, B1, C1, T2)
    Pvs_2 = antoine (A2, B2, C2, T2)
    Pa_cal = Pvs_1 * x_1 * gamma_1 + Pvs_2 * x_2 * gamma_2
    f = (Pa - Pa_cal) ** 2
    return np.sum(f) 

T2 = minimize(equilibrium, T2, method = 'L-BFGS-B').x

# 求y1
Pvs_1 = antoine (A1, B1, C1, T2)
tau12 = a[0] + a[1] / T2
tau21 = a[2] + a[3] / T2
gamma_1 = NRTL(x_1, x_2, tau12, tau21, alpha)
y_1 = Pvs_1 * x_1 * gamma_1 / Pa

plt.rc('font',family='Times New Roman') #将全局的字体改为“Times New Roman”的形式
plt.rcParams['xtick.direction'] = 'in' #将坐标轴设置为朝内
plt.rcParams['ytick.direction'] = 'in' #将坐标轴设置为朝内
plt.plot(x_1,T2,label='NRTL')
plt.plot(y_1,T2,label='NRTL,')
plt.plot(x_1_exp,T1,'^',label='exp')
plt.plot(y_1_exp,T1,'>',label='exp')
plt.xlabel('component 1 fraction')
plt.ylabel('T/K')
plt.legend()
plt.show()